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Abstract. Using a field-theoretic approach, we derive the first few coefficients of the exact low-density 
( "virial" ) expansion of a binary mixture of positively and negatively charged hard spheres (two-component 
hard-core plasma, TCPHC). Our calculations are nonperturbative with respect to the diameters d+ and 
d- and charge valences q+ and q~ of positive and negative ions. Consequently, our closed-form expressions 
for the coefficients of the free energy and activity can be used to treat dilute salt solutions, where typically 
d+ ~ d- and q+ ~ g_, as well as colloidal suspensions, where the difference in size and valence between 
macroions and counterions can be very large. We show how to map the TCPHC on a one-component hard- 
core plasma (OCPHC) in the colloidal limit of large size and valence ratio, in which case the counterions 
effectively form a neutralizing background. A sizable discrepancy with the standard OCPHC with uniform, 
rigid background is detected, which can be traced back to the fact that the counterions cannot penetrate 
the colloids. For the case of electrolyte solutions, we show how to obtain the cationic and anionic radii as 
independent parameters from experimental data for the activity coefficient. 



PACS. 82.70.-y Disperse systems 
electrolytes, molten salts, etc. 



52.25.-b Plasma properties - 61.20.Qg Structure of associated liquids: 



1 Introduction 

The two-component hard-core plasma (TCPHC) has been 
used for a long time as an idealized model for electrolyte 
solutions. In this model, also known as "primitive mod- 
el," the ions are spherical particles that interact with each 
other via the Coulomb potential and a hard-core potential, 
which avoids the collapse of oppositely charged particles 
onto each other. In the general asymmetric case, the pos- 
itive ions have a charge q+ e (where e is the elementary 
charge) and ionic diameter d+, while the negative ions 
have a charge — g_e and diameter d_. The particles are 
immersed in a structureless solvent whose presence is felt 
only through the value of the dielectric constant of the 
medium, and the system is (globally) electroneutral. 

The TCPHC in this formulation is a gross simplifica- 
tion of real systems. When two ions are at distances of 
the order of the size of the solvent molecules, the assump- 
tion of a continuum solvent breaks down, giving rise to 
so-called solvation forces [0. Such effects also lead to non- 
local contributions to the dielectric constant. Since the 
dielectric constant of ions or colloids typically differs from 
the surrounding aqueous medium, one also expects disper- 
sion forces to act between these particles. All these effects, 
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which contribute to "ion-specific" effects |||, and are pos- 
sibly more relevant than previously thought ||, arc not 
accounted for in the TCPHC, which treats the solvent — 
usually water — as a structureless medium within which 
the charged particles are embedded, neglecting the molec- 
ular arrangement that occurs around the ions. 

Nevertheless, even with such simplifications, the TCP- 
HC is far from being amenable to an exact treatment. A 
better understanding of this model is a necessary step if 
one wishes to develop more realistic approaches to charged 
systems. In this article, we turn our attention to the low- 
density, or virial, expansion of the TCPHC. Since we have 
exact results for the thermodynamic properties at low con- 
centrations using field-theoretic methods, we can obtain 
useful information on dilute systems with otherwise ar- 
bitrary hard-core radii and charge valences (like colloidal 
suspensions). In particular, we can test various approxi- 
mations to treat the TCPHC, like the mapping of a col- 
loidal solution on an effective one-component plasma. 

Due to the long-range character of the Coulomb po- 
tential, it is not easy to obtain the thermodynamic be- 
havior of the TCPHC through the usual methods of sta- 
tistical mechanics. For example, it can be shown ||,||,|{| 
that the straightforward application of the cluster expan- 
sion to the TCPHC leads to divergent virial coefficients. 
Mayer Q proposed a solution to this problem through an 
infinite resummation of the cluster diagrams, carried out 
such that the divergent contributions to the virial expan- 
sion are canceled. With this, he was able to obtain explic- 
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itly the first term in the virial (or low-density) expansion 
that goes beyond the ideal gas, which turns out to be the 
well-known Debye-Hiickel limiting lawQj. Haga[|| carried 
the expansion further and went up to order 5/2 in the 
ionic density. More or less at the same time, Edwards B 
also obtained the virial expansion of the TCPHC by mix- 
ing cluster expansion and field theory. In both cases only 
equally-sized ions were considered. 

The aforementioned methods typically depend on draw- 
ing, counting and recollecting the cluster diagrams which 
give finite contributions to the expansion up to the desired 
order in the density. This can be quite a formidable task, 
and unfortunately it is easy to "forget" diagrams that are 
relevant to the series (see for instance the comment on 
pp. 222-223 of Ref. ^|). Also, the generalization to ions 
with different sizes is rather complicated [||. Besides, the 
final results are typically not obtained in closed form, i.e., 
the final expressions depend on infinite sums that usu- 
ally have to be evaluated numerically, which reflects the 
infinite diagrammatic resummation. 

We generalize here a novel field-theoretic technique p0[, 
introduced for the symmetric TCPHC (q+ = q- and d+ = 
d-). We obtain the exact low-density expansion of the 
asymmetric (both in size and charge) TCPHC. This method 
does not use the cluster expansion (and resummation) and 
yields analytic, closed-form results. We go up to order 
5/2 in the volume fraction of a system where the sizes 
and the charge valences of positive and negative ions are 
unconstrained, that is, the results we obtain can be ap- 
plied, without modifications, to both electrolyte solutions 
(where anions and cations have approximately the same 
size and valence) and to colloidal suspensions (where the 
macro- and counterions have sizes and valences that can be 
different by orders of magnitude). In the colloidal limit, we 
demonstrate how to map the TCPHC on an effective one- 
component hard-core plasma (OCPHC) and obtain the 
corrections due to the exclusion of the background (i.e. the 
counterions) from the colloidal particles. For electrolytes, 
we obtain effective ionic sizes in solution, or thermody- 
namic diameters, from experimental data for the mean 
activity. The method developed here allows, in principle, 
to measure the diameters of cations and anions as indepen- 
dent variables, given that the range of experimental data 
extends to low enough densities such that the expansion 
used here is valid. It should be noted that the diameters 
obtained in this way are fundamentally different from the 
hydrodynamic diameters. 

This article is organized as follows. In section || we 
describe in detail the steps that lead to the low-density 
expansion. Readers which are not interested in this deriva- 
tion can go directly to section [| where we apply the ex- 
pressions obtained to two particular problems, viz., (i) col- 
loids, where one of the charged species is much smaller and 
much less charged than the other one and (ii) the mean ac- 
tivity coefficient (which is related to the exponential of the 
chemical potential) of electrolyte solutions, which is avail- 
able from experiments and can be compared to our results. 
Finally, section M contains some concluding remarks. 



2 The method 

We begin our calculation by assuming a system with A^_|_ 
positively charged particles with charge valence q + and 
diameter d + , and N- negatively charged particles with 
charge valence g_ and diameter d_ . Global electroneutral- 
ity of the system will be imposed at a later stage of the 
calculation. The canonical partition function Z is given 

by 



f 



N,\NJ. 



N+ dr (+) N - drH 



3=1 T '- 



where Xt.+ and A^._ are the thermal wavelengths, and 



rj + ' and ; are the positions of positively and negatively 
charged particles. The Hamiltonian TL is given by 



,(-) 



drdr' p Q (r)w ct/3 (r-r')p /3 (r') 



+ i| drdr' [ / 5 + (r)- / )_(r)] Vc (r-r')[ / 5+(r')-p_(r')] 

(2) 

where the charge-density operators of the ions are defined 
as 



Mr) 



M r ) = 9- 



£*(r-r< + >) 

i=l 
JV_ 

£<5(r-r«), 

i=l 



(3) 



and 6(r — r') is the Dirac delta function. The indices 
a and j3 in Eq. (3) stand for + and — , and the sum 
over a and /3 in Eq. (Q) runs over all possible permu- 
tations (viz. H — h, , — h and H — ), i.e., we consider a 

different short-ranged potential uj a p for each combination. 
The Coulomb potential is given by v c (r) — Ib/t, where 
Ib = e 2 / (Att e ks T) is the Bjerrum length, defined as the 
distance at which the electrostatic energy between two el- 
ementary charges equals the thermal energy ks T. Finally, 
-E S eif is the self-energy of the system 



E scl{ = ^±[uj ++ (0) + w e (0)] 



N-q 2 _ 



u;__(0) + « c (0)] 
(4) 

which cancels the diagonal terms in Eq. pi). 

We proceed by applying the Hubbard-Stratonovich trans- 
formation, the essence of which is given by 

e ~5 / drdr'p(r)t)(r-r')p(r') _ 

/D</>e _ 5 / drdr ' 0(r)f _1 (r-r')0(r')-«J dr0(r)p(r) 



JV(f)e~i / drdr ' 0(r)«" 1 (r-r')0(r') 



(5) 



where v(r) is some general potential and J T><j> denotes a 
path integral over the fluctuating field 4>. While this trans- 
formation can be used without problems when v (r) is the 
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Coulomb potential, for a short-ranged potential this can 
be more problematic: for instance, a hard-core potential 
does not even have a well-defined inverse function. We will 
anyway take this formal step for the short-ranged poten- 
tial, and, as we will see later, the way we handle the re- 
sulting expressions leads to finite (and consistent) results, 
viz., the virial coefficients 

Applying Eq. (||) to Eq. (|l|) we obtain the partition 
function in field-theoretic form 



makes the calculations much more difficult. Later, elec- 
troneutrality will be imposed order-by-order in the low- 
density expansion in a consistent way, any infinities aris- 
ing from the non-neutrality of the system will then be 
automatically canceled. 

In its full form, the grand-canonical partition function 
Q reads 



Z= [ V ^-^e-*°W + W- 



Q = 



(6) 



A- 



CXP (i^ J drh + (r)e 



-iq+4>(r) 



+ 



with the action 



dr/i_(r)e l<? -^W-7Y , (14) 



^ = \ E / drdr ' ^(r)^ - r')Mr') 

a,P=+,- J 

+ i J drdr' ^(r)^ 1 ^ - r')0(r'). (7) 

The inverse potentials are formally defined as the solution 
of the equation 

J2 f dr' u aP {v ~ r>^(r' - r ") = 8 a , 5{v - r") (8) 

13=+,- J 

(S a7 is the Kronecker delta) and 



J dr' v c (r - r')«jV - r ") = s ( r ~ r " 



(9) 



For the Coulomb potential, tt" 1 ^) = -V 2 5(r)/47r^ B . We 
also define 



1 



a «a[w aa (0)+« c (0)] /S 



dr 



-iq a [ip a {r)+a 0(r)] 



N„ 



and the normalization factors 

-| £ I drdr ' *Pa(r)u)~p{r-r' )ip{3{r' ) 



-i / drdr' 0(r)i; < r 1 (r-r')0(r') 



Z^p = J T>ip + Dip- e 
and 



(10) 
(11) 
(12) 



where TLq is given in Eq. (0). We defined the local non- 
linear operator 

2 

Mr) = exp(^[w QQ (0) +v c (0)] - ufcrfMr)) (15) 

(as before, a stands for both + and — ). At this point we 
rescale the fugacities such that A + = A + /\j, + and A_ = 
yl_/Ay _, i.e., the fugacities have from now on dimensions 
of inverse volume. 

Introducing the Debye-Huckel propagator, 

v^r - r') = ^(r - r') + I 2 S(r - r') (16) 

with the ionic strength 

h = q 2 + \+ + q 2 L\- (17) 

and after some algebraic manipulations, we finally ob- 
tain the grand-canonical free energy density. It is defined 
through g = — \n[Q)/V, and reads 

9 = -A + -A_-i/*(0)-i ta (^) 

_ ^ m ^ e A + /drQ + (r)+A_/drQ_(r)^ ^ 

where V is the volume of the system and the brackets 
(• • • ) denote averages over the fluctuating fields (f> and ip a 
with the propagators and 2dh is defined as 



--DH 



V<f)e 



-\ J drdr' (p(r)v D ^(r-r')^(r') 



(19) 



In order to make the calculations simpler we use the 
grand-canonical ensemble. This is achieved through the 
transformation 



and the functions Q+(r) and Q-(r) are defined by 



N + ,N-=0 



(13) 



where A + and A_ are, respectively, the fugacities (expo- 
nential of the chemical potential) of the positively and 
negatively charged particles. We perform the sum over N + 
and iV_ without constraints, i.e., without imposing the 
electroneutrality condition q+N + = g_iV_. Imposing this 
condition before going to the grand-canonical ensemble 



Q a (r) = h a (T)e-~ ^ (r) -l + ^ 2 (r)-i 9 > c (0). (20) 

In the preceding steps we obtained the exact expres- 
sion for the grand-canonical free energy density g, still 
without imposing electroneutrality. In order to obtain the 
low-density expansion of the free energy, we now (i) ex- 
pand g in powers of A + and A_ (up to a order 5/2), (ii) 
calculate the concentrations of positive and negative par- 
ticles and impose electroneutrality consistently, order-by- 
order, and (iii) make a Legendre transformation back to 
the canonical ensemble. 
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2.1 Expansion in powers of the fugacities 



and 



We start the expansion of g by noting that the Fourier 
transform of the Coulomb potential is v c (k) = Airls/h 2 . 
Using this, one is able to express v c (0) as 



«c(0) = 

It is easy to show that 



dk 4%£ B 
(2tt) 3 k 2 ' 



1 , Zdh 



dk 



In 1 



(2tt) 3 , s 
Using Eqs. @ and (§2|), one finds 



i/^ + Ilnf^S ] = J_ 
2 2 cW U I i 12tt 



4:ire B I 2 
k 2 



(21) 



(22) 



3/2 



Introducing the dimensionless quantity 



Av = \Utt£ 3 b I 2 



(23) 



(24) 



and expanding the last term in the rhs of Eq. ( |iq ) in cu- 
mulants of Q+(r) and Q-(r), one obtains the low-fugacity 
expansion 



X+Z 1 



x 2 

2 2 



^^--A + A_Z+-+0(A 3 ), (25) 



where the coefficients are given by 



Zi 



1 

r 



dr<Q+(r)) 



(26) 



and an analogous formula for Z 1 , and 
1 



Z^ = 



V 



drdr 



'{(Q + (r)Q + (r')>-(Q+(r)><Q + (rO>} 

(27) 

and similar formulas for Z 2 and ■ The symbol 0(A 3 ) 
in Eq. (|2"q ) means that any other contribution to the ex- 
pansion will be of order 3 or higher, i.e., with terms like 
X\, A+A, etc. Clearly, the expectation values in Eqs. ( p6|) 
and (g7|) contain additional dependencies on the fugacity 
A a via the DH propagator, Eq. (|l6|) , but, and this stands 
at the very core of our method, all expectation values can 
itself be expanded with respect to A Q and have finite val- 
ues as A Q — > 

In order to do a full expansion of Eq. ( p5| ) , one needs 
first to calculate the coefficients Z±, etc., in Eq. (p5[), for 
which the averages given in the appendix are needed, cf. 
Eq. (^| H|). Since v c (0) — «dh(0) = Avq, wc obtain 



zi 



1 



1 



(28) 



Zf = I drl e l<£W P ]**/3 



-la 9/3 [w a(3 (r)+a/3 udhM] _ ^ 



x 2 2 2 i \ 
^<7a<7/3«DH( r ) 



1 _ e 9o^o/ 2 _ e «|^o/ 2 



Note that 



l _B_ -Av r/ti 



(29) 



(30) 



was defined (through its inverse function) in Eq. (jlq). We 
now introduce the hard-core through the short-range po- 
tentials 



f+oc ifr^ + ^/2, 
otherwise, 



where the indices a and (3 stand again for + and — ; <i+ 
and d— are respectively the (effective) ionic diameters of 
the positive and negative particles. 

Since the expressions for Z^~ , etc., do depend on A+ 
and A_, and in order to have a consistent expansion of 
g in the fugacities, one has to expand Eqs. (|28|-|29|) in 
powers of A+ and A_ , up to the appropriate order, before 
inserting them into g, Eq. (^5|). By doing this consistently 
up to order 5/2 in the fugacities, one obtains the rescaled 
grand-canonical free energy density 



9 = d%g = 

— A + — A_ — miX 2 , — m 2 X_ — m^X + X^ — 



ni\ 2 + + 



n 2 X_ + n 3 A + A_ 



Av - 



ln Avq 



ri\ + +r 2 X 2 _ +r 3 A+A_ 



Avq In Avq — 



+ s 2 A_ 



M- 



-Av < Qi 



t + t 1 X++t 2 X- 
2q + q_ 



Avl + n 



q+x+-q-X- 



-n 



X+X- 



where we used dimensionless fugacities 



and 



A_[_ = cftX-\- 



X- = d±X-. 



qlXl + qtXl 
)+0(A 3 ), (32) 

(33) 
(34) 



In this expansion, we utilized that Avq scales like A 1 / 2 . 
The coefficients rru , etc. are given explicitly in the ap- 
pendix, Eqs. (HJ9|). 
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The coefficients and Q\ are given by the divergent 
integrals 



(35) 



These terms are present in Eq. (|32j) because global charge 
neutrality has not been yet demanded. By imposing this 
condition, these divergent terms cancel exactly, as is shown 
next. 



order by order. For instance, at linear order in A + , 
expanded form of Eq. (|3^) reads 



the 



14 



a <j- 



= 0, 



(40) 



f2 a = 2n£ B j drr 
o 

oo 

Ql = 2vr / drr 2 . 



which naturally gives ao = q+ /g_ . With the knowledge of 

ao, one can solve the next-order term (in this case X 3 / 2 ) 
and obtain the value of ai, and so on. The resulting coeffi- 
cients ao up to a$ are given in the appendix — cf. Eqs. (|93| 

As this order-by-order neutrality condition is imposed, 
one notices that the terms J?o and fi\ in Eq. (|3^) are ex- 
actly canceled in a natural way, without any further as- 
sumptions. The resulting expression for g, now expanded 
only in one of the fugacities (in this case A+), is then a 
well behaved expansion (we omit this expression here since 
it is quite lengthy). With this, one can finally obtain the 
canonical free energy (as a density expansion) through a 
Legendre transform. 



2.2 Imposing electroneutrality 

The electroneutrality condition ensures that the global 
charge of the system is zero, i.e., q+N + = q_N_. In the 
grand-canonical ensemble, N + and iV_ are no longer fixed 2.3 The canonical free energy 
numbers but average values. This means that the elec- 
troneutrality condition in the grand-canonical ensemble is 
given by 

q+(N + ) =q-(N_). 



Defining the rescaled ion density 



c+ = dl(N + )/V, 



it is easy to show that 



c+ = — A-| 



<9Aj 



(36) 



(37) 



(38) 



In order to obtain the free energy in the canonical ensem- 
ble, we use the back-Legendre-transformation 

f = g+Z + ln(x+) +£_ln(A_), (41) 

where the dimensionless canonical free energy density is 

f = d\F/Vk B T (42) 

(F is the canonical free energy). Note that at this point 
A_ is a function of A + , according to Eq. (|39|). 

The first step to obtain / is to invert the expression 



with an analogous formula for c_ ; c+ is the volume frac- 
tions of the positive ions[^2|. As one imposes Eq. (|36|), the 
fugacities will depend on each other in a non-trivial way 
such that the system is, on average, neutral. If the system 
were totally symmetric (i.e., q + = a_ and d + = d-), this 
dependence would be given by the relation A + = A_ ||To||. 
However, this is not the case here: one needs to find the 
relation between the fugacities order-by-order. 

First, assume that A_ can be expanded in terms of A+ 
such that 



given in Eq. (38) such that A+ is obtained as an expansion 
in c + , neglecting any terms of order cr^ or higher. With 

this, we obtain A+ = A + (c + ): plugging this into Eq. (fill) , 

we finally obtain /, which reads 



/ - fid 



, R ^/ 2 



B 2 i og ci In c+ 



B 5 /2\o e c+ 2 lnc+ 



0(cl) 



~ ~ — 3/2 — 2 ~2 ~ 

A_ = aoA + + a±X i + a 2 X + + a^X + In A+ 



+ <24 Ai/ 2 + 0,5 A V 2 In A^ 



+ 0(A 3 + ). (39) 



This is inspired by the expanded form of the grand-canonical 
free energy g, Eq.@. ^ 

After calculating c + and c_ from g, Eq. (|32|), according 
to Eq. (|38|) , and insertion into the electroneutrality condi- 
tion Eq. (|3|), we substitute the fugacity A_ by its ex- 
panded form Eq. (|39|). This leads to the expanded form 

(up to order X 5 / 2 In A + ) of the electroneutrality condition. 
Solving it consistently, order-by-order, yields the values of 
the coefficients in Eq.(|39|), which ensure electroneutrality 



Defining the valence ratio parameter 

V = <7-/q+, 
the diameter ratio parameter 

C = d-/d+, 
and the coupling parameter 

e+ = q 2 + £ B /d+ 



(43) 



(44) 



(45) 



(46) 



(which is the ratio between the Coulomb energy at contact 
between two positive ions and the thermal energy k B T), 
the coefficients in Eq. (J43l) can be explicitly written as 



/id = c + In c + + — lnf — 



1 



(47) 
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Fig. 1. Coefficients (a) Bi and (b) B5/2 of the free energy expansion Eq. ([43|) as functions of e+ for different values of r] = q-/q+ 
(with d+ = d- or £ = 1). Notice that B2 is normalized by the second virial coefficient -B^ of a pure two-component hard-core 
gas. 



which is the ideal contribution to the free energy, 



Sdh = -|\/7re^[l + ry] 3 



(48) 



which is the coefficient of the Debye-Hiickel limiting law 
term (order 3/2 in c+). It is useful to define the function 

H(x) = —-2-f + -^e- x \2-x + x 2 ]-r(0,x)-\nx. (49) 
6 x A 

where 7 is the Euler's constant and r(a,b) is the incom- 
plete Gamma- function |i~3| . Using H(x), the higher order 
coefficients can be explicitly written as 



B 2 = --eU-H(e + )-\ne + + 2 V 2 



2 V e A 



V l + eJ 



In 



H 



rj 2 e+ 



In 



7? 2 e H 



2t7 2 [l - r/ 2 ] In 77 + - [l - r) 2 ] ln(367re 3 [1 + 



(50) 
(51) 



%2 = § K] ^ [1 + V] 1/2 1 f + H(e + ) + ln(e + ) 



In 



5 + Fff 



■ [1 + 77] 2 [5 - 1277 + 17r? 2 - 12?7 3 + 5?7 4 ] 



2t7 3 [1 + 77 3 ] In 77- i [1 + r? 3 ] 2 lnfWe* [l + 77] ) 



(52) 



and 



5/21og 



lK] ,/, [i+i] 1,J [i+i'] i 



1 1/2 r 



(53) 



The free energy Eq. ( Eg ) is the exact low density ex- 
pansion of the asymmetric TCPHC and constitutes the 
main result of this paper. The only parameter that is de- 
manded to be small is the ion density c+, that means, this 
result is non-perturbative in the coupling e+, charge ratio 
77 and size ratio £. 

We chose the positive ions as the "reference species" 
(i.e., the expansion is done with respect to c+) without any 
loss of generality, since the relation between c+ and c_ is 
fixed through the electroneutrality condition. As consis- 
tency checks, we notice that our expression for / is sym- 
metric, as expected, with respect to the simultaneous ex- 
change of d+ with d- and q + with q_. Also, in the limit 
d + = <i_ and q + = q_, we recover the same expression 
as previously calculated in Ref. M] for totally symmetric 
systems. Finally, as one turns off the charges in the sys- 
tem (or equivalently, as one takes the limit e + — > 0), the 

pure hard-core fluid is recovered, i.e., / becomes the usual 
virial expansion with -B3/2 and B 5 / 2 equal to zero and B 2 
given by the second virial coefficient of a two-component 
hard-core gas, B\^ c = B 2 (e + — > 0), which in our units 
reads 



2£ 3 , [1 + £] S 



277 



(54) 



This limit can be also understood as the high-temperature 
regime: as the thermal energy largely exceeds the Coulomb 
energy at contact, the hard core interaction becomes the 
only relevant interaction between the particles. 
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Fig. 2. Coefficients (a) B2 and (b) -B5/2 of the free energy Eq. ( f43| ) as functions of e + for different values of £ = (with 
g+ = g_). Notice that £?2 is normalized by the second virial coefficient B^ c of a pure two-component hard-core gas. The inset 
in (a) shows the behavior of B2 close to e+ = on a different scale. 



3 Results 

3.1 The virial coefficients 

The behavior of the coefficients B2 and J5 5 / 2 as functions 
of the coupling parameter e + are depicted in Figs. [j]and |], 
the behavior of i?2iog and B 5 /2i os is rather trivial and not 
shown graphically. In Fig. ^ the ionic diameters of positive 
and negative ions are kept equal, d + — d_, and the ratio 
between the charge valences, ry, is varied, while in Fig. || 
the charge valences are equal and the ratio between the 
ionic diameters is varied. These figures highlight the fact 
that both coefficients diverge as e+ goes to infinity. In this 
limit, 

7T[l+£] 4 1 / 277£_ 

T exp 

4 if 6+ 



Bo 



iff!) 



(55) 



and 



B 



5/2 



7T 3 / 2 [i + ^]4yr+^ 



' exp lTTi> 



(56) 



Note that B 5 / 2 diverges faster (and with opposite sign) 
than B2 , reflecting that this low-density expansion is badly 
converging, as we will discuss further below. The exponen- 
tially divergent behavior of B2 and -B5/2 when e+ — > 00 
is due to the increasing importance of the interaction be- 
tween oppositely charged particles (ionic pairing) as the 
coupling parameter increases JIJ,^5|,|l^], corresponding for 
instance to lower temperatures. This is confirmed by not- 
ing that the argument in the exponential occurring in both 
asymptotic forms Eqs. ( |55| ) and (M), viz. 2 77 e + /[l +£], 
can be re-expressed as 2 q + q_ £B/\d + + <2_], which is the 
coupling between positive and negative ions (in this case, 
the ratio between the Coulomb contact energy between 
oppositely charged ions and the thermal energy fcgT). 

To estimate roughly the ionic density c+ up to which 
the expansion in Eq. ( j43| ) is expected to be valid, we use 
the following simple criterion: the terms proportional to 
-Bdh and -B5/2 give the same contribution to the free en- 
ergy when |-Bdh|<x|/ 2 = | -S5 /2 1 c^/ 2 or, in other words, at 



a critical density c+ = \B-£>n/B 5 / 2 \- Note that we set up 
this criterion separately for the integer and fractional co- 
efficients, since the scaling behavior for the two different 
classes is very different and no meaningful results can be 
obtained by mixing them. The analogous criterion for the 
integer terms leads to a critical density c+ = {l+X/ri)/\Bo\ 
where for the ideal contribution in Eq. |47| we replace the 
logarithmic term by a linear one (which is the ideal con- 
tribution to the osmotic pressure) . We tacitly assume that 
the higher-order terms (which we have not calculated) 
show the same relative behavior, an assumption which 
seems plausible to us but which we cannot check. In Fig. 
H we show the estimate for the critical densities obtained 
from a) the ratio of the integer terms and from b) the ratio 
of the fractional terms. In a) the critical density settles at 
a finite value for vanishing coupling parameter e + — * 0, 
and decreases to zero for increasing coupling parameter 
(the divergence in the curves at finite value of e+ is not 
significant since it is caused by a change in sign of i?2-)- 
This means that for small coupling constants the integer 
coefficients are expected to show regular convergence be- 
havior for not too large concentrations. For large values of 
c + on the other hand, the range of densities which can be 
described by the expansion quickly goes to zero. In b) the 
displayed behavior is more complex. The expansion of the 
coefficient B 5 / 2 for small values of the coupling parameter 
e + leads to 



B 



5/2 



= ^/ 2 yir^(e~i) 2 (e + i> 



.3/2 



(57) 



plus terms with scale as higher powers in e + . The impor- 
tant point is that except for the singular case £ = 1, that 
is for particles of identical radius, the leading term of -B5/2 
3/2 

scales as e + and thus identical to the lower leading term 

5 /2 

£?dh- In the case £ = 1 the leading term scales as e + 
and thus shows a different scaling behavior than -Bdh- In 
Fig. ||b) this difference is cleary seen: For £ = 1 the curves 
diverge as e + — > and the convergence of the fractional 
density terms in the series is expected to be guaranteed. 
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Fig. 3. Ratios of the coefficents in the low-density expansion, which serve as rough estimates of the maximal rescaled ion 
density £+ up to which the expansion is valid (see text). 



For £ = 0.1 the behavior is similar to the results in a), 
showing a saturation at a finite value as e+ — > (and 
an unsignificant divergence at finite £+). All curves go to 
zero as the coupling parameter grows. In conclusion, the 
reliability of the low-density expansion becomes worse as 
the coupling parameter increases, but one can always find 
a window of small densities within which the expansion 
should work. 

Finally, we mention that in a different calculation scheme, 
based on an integral-equation-procedure within the so- 
called mean-spherical-model (MSM) approximation, sim- 
ilar results for the free energy and other thermodynamic 
functions have been obtained p7| , [l^ ,pj|. Those results also 
reproduce the limiting laws, namely the hard-core behav- 
ior as the coupling parameter goes to zero, and the leading 
Debye-Hiickel correction at low density, and give a quite 
satisfactory description even for much larger values of the 
density. It is important to note, however, that the next- 
leading (beyond Debye-Hiickel) terms in the low-density 
expansion are not correctly reproduced by the MSM ap- 
proximation. 



With this in mind, let us take the following limit: as- 
sume the parameters describing the positive ions (repre- 
senting the colloids) d +1 q + and c+ fixed, and make both 
rj = (?-/<7+ and £ = d-/d + vanishingly small (this has to 
be done with some care, since the limit rj = is not well- 
defined) . One can then rewrite the free energy Eq. ( [£31) up 
to second order in c+ as 



3.2 The colloidal limit 

In colloidal suspensions, flocculation or coagulation (driven 
by attractive van der Waals interaction between colloidal 
particles) can be prevented by the presence of repulsive 
electrostatic forces. These suspensions are typically quite 
dilute, with volume fractions of colloidal particles usu- 
ally not higher than a few percent. The macro-particles 
normally have dimensions pOl Eh, E2| ranging from 10 to 
1000 nm and charges of several thousands e (elementary 
charge) , with much smaller counterions that have a charge 
of a few e. In such systems the charge and size asymme- 
try between ions and counterions is immense; the TCPHC 
with unconstrained charge and size asymmetry is a suit- 
able model for dilute colloidal suspensions, and the free 
energy Eq. ([l3|) can be used to study the thermodynamic 
behavior of such systems. 



r = /id 



-B 



J -3/2 
DH C + 



(B* G ' cl + B\ 



B: 



;1 

21o g c + 



log C 4 



o{tl 2 ) 



where /id is the ideal term Eq. (p7|) and 



B, 



HC.cl 



2e t [i+£] 3 



2r) 



(58) 



(59) 



is the hard core contribution due to the counterion-counterion 
and macroion-counterion interactions only, without the 
macroion-macroion contributions, which explains the dif- 
ference to the full hard-core virial coefficient in Eq.(|54|). 
It is necessary to treat this term separately, as it does 
not have a well-defined behavior in the double limit rj — > 
and £ — > 0. The other coefficients in Eq.([58|) follow 
from Eqs.(||), @), (§!]) by performing the limits Bg H = 
Bdh(?? - CU - 0), Bf = B 2 ( V - 0,£ -> 0) - B^ c \ 



DC! 



= -B2io g (?7 — > 0, £ — > 0), and are given by 



R cl _ 



+ ■ 



R cl _ r c 3 



{ff(e+) -iln(36^+)} + 



r e +» 



Del 

-°2W 



6 



(60) 



(61) 



(62) 



HC,cl 



It is to be noted that the contributions in and B. 
contain terms that scale with 1/rj. For dilute colloidal sys- 
tems therefore the ideal contribution of the counterions 
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Fig. 4. Coefficients Bf (Eq. (|65|)) and d 2 (pure OCPHC from 
Ref. JTTJ) as functions of the coupling e+. Notice that both 
coefficients are normalized by the second virial of a pure one- 
component hard-core gas. 



dominates the free energy and cannot be neglected. Ef- 
fects due to the electrostatic interaction between the ions 
are corrections to the ideal behavior. 

At this point we briefly introduce yet another model 
that is also widely used to describe charged systems. It 
is the one-component plasma (OCP), which in its sim- 
plest form consists of a collection of N equally charged 
point-like particles immersed in a neutralizing background 
that ensures the global charge neutrality of the system 
(in the TCPHC electroneutrality is ensured by oppositely 
charged particles). The OCP, or its quantum mechanical 
counter-part ("jellium") has been used in different con- 
texts in physics, as for instance to describe degenerate 
stellar matter (the interior of white dwarfs or the outer 
layer of neutron stars) and the interior of massive plan- 
ets like Jupiter. Another example comes from condensed 
matter physics, where jellium is often used as a reference 
state when calculating the electronic structure of solids. 
For reviews see Refs. jgj| |4[gg] . 

When the particles have a hard core, the OCP is called 
one-component hard-core plasma (OCPHC): what we will 
see next is that the electrostatic contribution to the free 
energy in dilute colloidal suspensions can be almost de- 
scribed through the OCPHC. If we compare the coeffi- 
cients Eqs. (|60|~|5l|) with the ones previously obtained JlO[ 
for the low density expansion of the OCPHC 

/ocphc = c In c— c+ d 3/2 c 3 / 2 + d 2 ? + d ln 2 7? In c+ 0(e> /2 ) , 

(63) 

where c is the volume fraction of particles in the OCPHC, 
we find (cf. Eqs. (51), (56) and (57) of Ref. @) 



^DH — ^3/2 ! 



(64) 



d 2 + -e H 



and 



pel 



(66) 

where the notation used in Ref. |l(| is kept in Eq. (|63| ) and 
on the rhs of Eqs. (]64]-|66|). The comparison between Bf 



and d 2 is shown in Fig. where we rescaled both coeffi- 
cients by their value at vanishing coupling, d 2 {e+ = 0) = 
B^(e + = 0) = 2ir/3. This shows that, including the or- 
der lnc+ (i.e., for very dilute colloidal suspension), the 
OCPHC is almost recovered, except an additional term in 
Eq.(j35j). This additional term can be quite important at 
intermediate values of e+ since it changes the sign of the 
virial coefficient, as seen in Fig. 3. 

This is in fact what one would intuitively expect: each 
macroion has around it a very large number of small neu- 
tralizing counterions which act like a background. Our re- 
sults show that the background formed by counterions, 
which is a deformable background since the counterion 
distribution is not uniform, acts almost like a rigid homo- 
geneous background, as it is assumed in OCP calculations. 
However, there is a small difference between the TCPHC 
in this limit and the OCPHC: in the latter, the background 
penetrates the particles, while in the TCPHC it cannot 
(for a thorough discussion of this difference see [p6[). In 
our calculation, this is reflected in the term, where 
7re + /2 is the positive extra cost in the free energy that the 
OCPHC has to pay (at this order) to expel the background 
from the hard-core particles. The first consequence of the 
non-penetrating background is that the effective charge of 
the colloids increases by exactly the amount of background 
charge that is expelled from the colloidal interior. A sim- 
ple calculation shows that the increased effective charge 

of colloids for the non-penetrating background turns 
out to be 

„eff _ Q+ 



1 — 7TC+ /6 



(67) 



In the low-density expansion of the OCPHC, Eq. (p3|), one 
would have to reexpand all coefficients with respect to the 
density, using Eq.(|37|). However, since the leading term 

3 /2 

which depends on the charge q + goes like c + , the ef- 
fect would come in at order c+ 2 and therefore is not 
responsible for the additional term in Eq.(|65|). The rea- 
son for the extra term in Eq.(|65|) has to do with an in- 
crease of the self-free-energy of a colloidal particle, which 
can be understood in the following way: Assume that the 
OCPHC is very dilute, such that each colloidal particle 
and its neutralizing background form a neutral entity (in 
the spirit of the cell model, see for instance Ref. |27|) that 
can be regarded independent from the other particles. The 
free energy difference per colloidal particle between a sys- 
tem without penetrating background and with penetrat- 
ing background (Fig. |5j) is then given by two contribu- 
tions: one coming from the entropy lost by the background 
(formed by AL. counterions) since it cannot penetrate the 
macroions, given by 



AF ei 



(65) N+k B T 



-^ln(l-7r£ + /6) = ^^± + 0(4). (68) 



This is one of the contributions to the second virial co- 
efficient due to the hard-core interaction, and already in- 
cluded in Eq. (59). This term does not depend on the 
coupling e +) and therefore has nothing to do with the ex- 
tra term in Eq.(|65|). The second contribution to the self- 
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Fig. 5. A colloidal particle (with diameter d+ and charge va- 
lence q+) in the OCPHC model (a) with and (b) without the 
penetrating background. 



energy is electrostatic in origin. Defining the charge distri- 
butions for the situations in Fig. 4 as p a (r) = q+S(r) — p 
when the background can enter the colloids and p&(r) = 
q+6(r) — po [l — 8(d + /2 — r)] when the background cannot 
enter the colloids (where po = c + q + is the background 
charge density and 6(x) = 1 if x > and otherwise), the 
electrostatic self-energy difference between the two cases 
reads 

_AF, 

N. 



~^tr = \\ drdr ' { Pb ^ Vc ( r ~ r ')pe>( r ') 

- p (r)u c (r - r')p a (r') 



poq+ATrlB 



d+/2 



dr r 

7T 



+ 0(c 2 + ) 

c+e+ + 0{c 2 + ). (69) 



It follows that the free energy difference per volume is 
given by 

c + AF el 7rcle+ 



(70) 



where A+ and A_ are respectively the fugacities of the pos- 
itive and the negative ions. The mean activity coefficient 
f± is the ratio between the mean activity of the electrolyte 
and that of an ideal gas, in general given by 



f± = exp 



dfe 



q+ + q- dc+ J ' 



(72) 



N + k B T 2 

Note that this is a positive contribution to the c 2 term 
with the same coefficient as the extra term in Eq. (p5|). 
Therefore, we conclude that this extra term is due to the 
electrostatic energy associated with expelling the back- 
ground from the colloidal volume. This is in accord with 
more general results on the difference between OCPHC 
models with penetrating and excluded neutralising backgrou 
In summary, whenever using the OCPHC to describe highly 
asymmetric charged systems such as colloids, one has to 
to take into account the exclusion of the background from 
the macroions. As we demonstrated, if this is taken into 
account, then the TCPHC maps exactly (at least up to 
order ci) onto the OCPHC. 



3.3 Ionic activity and diameters 

In electrochemistry, it is usual to define the mean activity 
A± of an electrolyte as 



A j 



Al-A 9+ 



!/[?++?-] 



(71) 



for a two-component system where the positively charged 
particles are used as reference species in the same way as 
in Eq. ([43|). / ex , the excess free energy, is the difference 
between the free energy of the interacting systems and 
the free energy of an ideal gas, i.e., / ox = / — fid- 

There are different ways of measuring f±, as for in- 
stance, through the change of the freezing point of the 
solvent (usually water) with the addition of saltp8[, by 
measuring the change of the potential difference between 
the electrodes of a concentration cell as salt is added 
psf (potentiometry) , or by direct measurement of the sol- 
vent activity through vapor exchange between a solution 
with known activity and the sample p0||2^ ] (isopiestic). Al- 
though dating from the early nineteen hundreds, these are 
still the most common techniques used today, especially 
the potentiometry, which is regarded as the most precise 
technique of all. The values of f± are tabled as function of 
the salt concentration for many different electrolytes [|3l"[ 

HID- 

From the free energy Eq. (^3|) and the definition Eq. ( [72] ) 
we can obtain the low density expansion for the mean ac- 
tivity coefficient of a q + : g_ salt. To compare with ex- 
perimental results, it is useful to note that c+ = 6.022 x 
IQ^q-dSg, where d + is in Angstroms and g is the salt 
concentration in moles/liter. After the appropriate expan- 
sion, f± reads 

f ± = 1 + ^dh£> 1//2 + v\q + vi\ og g\ng + ^ 3 / 2 £> 3/2 

+ V3/2io g Q 3/2 li-iQ + 0(g 2 ) (73) 

(the order 3/2 in the mean activity coefficient is the one 
consistent with a free energy given up to 5/2). The coef- 
ficients ^dh, v li etc are given in the appendix, Eqs. (|99|- 
103). Experimentally determined activity coefficients are 
typically measured at constant pressure, while the the- 
oretical results are obtained for constant volume. Note 
fclgp^i the limit of vanishing density, the distinction be- 
tween the activities calculated in the Lewis-Randall or in 
the McMillan-Mayer descriptions (constant pressure ver- 
sus constant volume ensemble) becomes negligible in 
comparison to the first corrections to ideal behavior. The 
experimental data at low density can be directly compared 
with the theory without the need to convert between the 
two ensembles. 

At infinite dilution, the mean activity coefficient Eq. (|7q) 
goes to 1, which is the prediction for an ideal gas. The first 
correction to the ideal behavior is the term vdhQ 1 ^ 2 , which 
is the prediction one obtains from the Debye-Hlickel lim- 
iting law (DHLL), and is independent of the ionic diame- 
ters, see Eq. (|9|). This means that there is always a range 
of concentrations where different salts (but with the same 
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this will determine uniquely the single unknown parame- 



Fig. 6. Experimental[||, |l| 

and fitted theoretical mean activ- 
ity coefficient f± for various salts as a function of the density 
q (in mole/liter), for various 1 : 1 salts, viz. HC1 (triangles 
and dot-dashed line, d — 4.0 A), NaCl (stars and dotted line, 
d — 3.7 A) and KC1 (diamonds and dashed line, d = 3.6 A). 
The full line denotes the Debye-Hiickel limiting law, which 
scales with g 1 ^ 2 . The fit was done using the data points shown 
(up to g — 0.02 mole/liter). 



q+ and q-) will deviate from ideality, but have the same 
activity. At higher concentrations, other terms have to be 
taken into account, and the ionic sizes begin to play an 
important role. In fact, as one fixes £ B 7.1 A in water 
at 25 °C) and the valences q+ and q_, the only free pa- 
rameters in the rhs of Eq. (|73|) are the ionic diameters d+ 
and d—, which can be used in the theoretical predictions 
to fit the experimental values. This leads to effective equi- 
librium values of the ionic diameters (when in solution), 
which we call the "thermodynamic diameter" , in contrast 
to the bare diameter pi (obtained through crystallographic 
methods) and the hydrodynamic diameter (obtained from 
mobility measurements |34| ) . By construction, it is this di- 
ameter which should be used in equilibrium situations if 
one wants to describe an electrolyte solution as a TCPHC, 
as for example computer simulations of electrolytes. 



3.3.1 Fitting assuming one mean diameter 

We now show the fitting procedure assuming that d + = 
d_ = d, where d is the mean diameter. We also restrict 
here this procedure to 1 : 1 salts, where more experimental 
data is available at reasonably low densities (below ~ 0.05 
mole/liter). For asymmetric salts, this method demands 
precise data at the range below ~ 0.01 mole/liter, where 
Debye-Hiickel is often assumed to account for all effects 
and few experimental points are available. 

The assumption of equal ionic sizes has been often 
used in the past to fit activities to modified Debye-Hiickel 
theories p8[ to account for the ionic sizes (as previously 
mentioned, the Debye-Hiickel limiting law is insensitive 
to it). Since we force the two diameters to be equal, we 
only need the expansion for f± up to linear order since 



ter. 



_5|,[n| we subtract 
the diameter-independent term and obtain the difference 



From the experimental data f ± 



Ai+ = £ 



_ rexp 



1/2 _ exp 



(74) 



The fit to this function leads to the coefficient of the linear 
term v^ xp ', and that can be used to determine the diameter 
by solving the equation Vi(d) = u± (with d as unknown). 
In Fig. || we show the experimental p5|, [3l| f± for HC1, 
NaCl and KC1 and the theoretical results (up to linear or- 
der in g) using, respectively, the diameters 4.0 A, 3.7 A and 
3.6 A. These values come from the fitting described above 
applied to the experimental data in the range g = 0.01 to 
0.05 mole/liter. The diameters obtained are very close to 
the ones obtained in Ref. ]35"| with a similar fitting, but 
using the DH theory with an approximate way to incor- 
porate the ionic sizes. 

The fact that cZhci > G?NaCi > c^KCi: which is the op- 
posite to the sequence of bare diameters, is a consequence 
of the existence of hydration shells around the ions which 
tend to be larger for smaller bare ion sizeQ]. Notice how- 
ever that the values obtained here lie between the bare 
diameters and the hydrated values available in the liter- 
ature. This is not surprising since the effective hard-core 
size simultaneously reflects both the presence of the hy- 
dration shells and the "softness" of the outer water layer 
as two oppositely charged ions approach each other. 



3.3.2 Fitting assuming two diameters 

If we now assume that both d + and e?_ are unconstrained, 
we have to use one more term in the expansion of the 
activity coefficient and solve the coupled equations 



exp 



Z/i(<i + , gL) = V 

"3/2(d + ,d_) 



exp 
= "3/2 



(75) 



where the coefficients v\ p and are obtained by fitting 
the function 



Ai± = i± p - v-ong 112 - ^Hogging 

= v\ xp g + v e 3 xp g 3/2 + 0(g 3 / 2 In g). (76) 

For the actual fitting we divide by the density and obtain 



exp 1/2 
y 3/2 £ • 



(77) 



Note that we subtracted the term proportional to g In g, 
since like the DH term it is independent of the ionic diam- 
eters. What is interesting about this fitting procedure is 
that the two ionic diameters are independent parameters, 
that is, the effective sizes obtained through this fitting do 
not depend on the size of a "reference" ion, contrary to 
what happens when using crystallographic methods p6[. 
The method we show here can in principle lead to use- 
ful results, as long as the experimental data is accurate 
enough at very low densities, as we now demonstrate. 
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Fig. 7. Mapping ionic diameters of a 1 : 1 salt in water at 
room temperature into coefficients of the rescaled activity coef- 
ficient Ai± defined in Eq. (jr^) . The inset shows different paths 
in the (d-;d+)-space (with the ionic diameters between 3 A 
and 27 A) which are m appe d into line s in the (fi;i/3/ 2 )-space, 
as determined by Eqs. (100) and (102). By knowing the values 
of v\ and ^3/2 f° r a certain salt at very low densities, one can 
make the "inverse-mapping" and extract the values of d- and 
d+. The solid line, the symmetric case with d+ = d_, is the 
lower envelope off all other lines which are obtained by varying 
one diameter while keeping the other fixed. 



Fig. shows the mapping between the parameter space 
of ionic diameters and the two coefficients v\ and ^3/2 (for 
a 1 : 1 salt in water at room temperature). The map- 
ping is restricted to the range 3 A < d + < 27 A and 
3A<gL<27A. This allows the "inverse mapping" be- 
tween the (^i;j/3/2) _ space and the (<i_ ;c?+)-space: any sys- 
tem with ionic diameters within the values above should 
have the coefficients v\ and ^3/2 inside the region in the 
(z^i;i'3/2)-space, as shown in Fig. ^. One problem becomes 
obvious: the coefficient ^ 3 / 2 is larger than the coefficient 
of the linear term, i>i, meaning that experimental data 
at very low concentrations are needed. As an example, 
Fig. I shows Af±/g for NaCl|l| as a function of g 1 / 2 , 
which should be asymptotically linear in the limit g — > 
and therefore give the coefficients v\ and ^ 3 / 2 by a simple 
linear fit according to Eq.(f77j). The lines shown are two 
possible asymptotic linear fits to the experimental points 
(line a is given by Af±/g = 2.59 — 9.94g 1 / 2 , and line b 
by 2.47 - 6.00£ 1/2 ). The inset to Fig. | gives the corre- 
sponding positions in the (ivi;i/ 3 / 2 )-space (cf. Fig.|^) for 
the two lines. Clearly, only one of the fits (line b) leads to 
reasonable values for the diameters (between 3 and 9 A as 
seen from the position of point b in the inset, d+ = 3.8 A 
and d- = 5.4 A as obtained by solving Eq. fl75|)), while 
the other fit (line a) leads to unreasonable values for the 
diameters (outside the range 3 to 27 A). In other words, 
the asymptotic extrapolation of the experimental data is 
very sensitive to small errors, and in order to obtain the 
diameters with this method one needs more accurate data 
at very low densities (which to the best of our knowledge 
is not available in the literature) . Although NaCl was used 
as example, the situation is identical for other salts. 



Af + 




Fig. 8. Fitting procedure (NaCl in water at room temperature) 
for two unconstrained diameters. Lines a and b correspond to 
two possible asymptotic fits of Ai±j g as g — > (see text). 
From these lines one can extract v\ and ^3/2, which allows 
determination of the ionic diameters (inset, which shows a part 
of the mapping done in Fig. [?]) Notice that while the two lines 
are reasonable asymptotic fits to the data, one (line b) leads 
to acceptable values for the ionic diameters, d+ = 3.8 A and 
cL = 5.4 A, while the other (line a) leads to values that lie 
outside the range 3 to 27 A (which is not reasonable). This 
means that more precise data in the range g < 0.01 mole/liter 
is needed to obtain the correct values for the diameters with 
this method. Experimental data from Ref. fell]. 



4 Conclusions 

Using field theoretic methods we obtained the exact low 
density ("virial") expansion of the TCPHC up to order 
5/2 in density. In its general form, the model can be ap- 
plied to both electrolyte solutions and dilute colloidal sus- 
pensions (when the van der Waals forces are unimpor- 
tant); the free energy derived here provides a unified way 
for handling both systems in the limit of low concentra- 
tion. As the calculations show, the generalization to short- 
ranged potentials other than hard core is possible. 

The behavior of the coefficients B2 and -B5/2 suggests 
that the scries is badly convergent, meaning that the in- 
clusion of higher order terms does not necessarily extend 
the validity of the free energy to larger densities. We saw 
that the divergent behavior of these coefficients is related 
to the ionic pairing[Q[l^,[l6| which is favored as the cou- 
pling increases (this is not present e.g. in the OCPHC@). 
This also means that such a low-density expansion is not 
very useful to study the phase behavior [|37| of ionic sys- 
tems: in the situation of phase-separation between a very 
dilute phase and a dense phase, typically the average den- 
sity in the dense phase already falls outside the range of 
validity of the low density expansion. 

In applying our results to colloidal systems we con- 
cluded that at low density the counterion entropic contri- 
bution dominates over the electrostatic contribution due 
to the macroions; the latter contribution can be described, 
in this limit and up to this order, by an OCPHC corrected 
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by effects due to the exclusion of counterions from the col- 
loidal particles. 

Finally, we used the theoretical results for the mean 
activity coefficient to fit experimental data and extract 
effective ionic sizes. In the simplest fitting, where we as- 
sumed the two ionic diameters to be equal, we obtained 
sizes that are reasonable and close to what one would ex- 
pect from the results obtained by other methods. For the 
more interesting case where the two ionic sizes are taken as 
free parameters and determined independently, we showed 
that one would need more experimental data for the mean 
activity coefficient at very low densities (which, to the best 
of our knowledge, is not available in the literature) to ob- 
tain the correct values for the diameters. With the proper 
experimental data it would then be a simple matter to ob- 
tain the effective thermodynamic ionic sizes, which could 
serve as useful input for computer simulations of two- 
component-plasmas with hard-cores to model electrolyte 
solutions. 
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Appendix 

Averages needed for Z\ and Z 2 



The following expressions were used to obtain Eqs. (28 



|29| ). In order to have more compact formulas we use the 
Greek letters a and /3 instead of + or — . For instance, 
aq a means both +q + and — q_. 



(MO) 



= e 9>=(0)/2 



-«>dh(0)/2 



(78) 



(79) 



m 2 = 77 



3 d% 



2n q%qti%{ „(2q+q_£ B 
m 3 = — — ^ — <! H ' 



d 3 



riq+q-Z B \ . / 3[d+ + d-] \ 
\d,+d_) n V 2l n ) 



3 d 3 



(^)-( 1 f) 



n 2 = 77 



T<£4f 5 n „f<£.l B 



3 dl 



2H 



In 



"3 



2nq 3 + q 3 L£ B ( 5 



79+9- 
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12d 2 _ 

B 



H 



(85) 
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(87) 
2q+q-£ B 



q+q- In 



2\d 4 



q\ + q 2 _ 



In 



d + + gL 
3[d+ + d. 
2£ R 



Pi = tt- 



:-! >r[ P2 3 d\ 

2tt q%£ 3 B 



7rq 6 .£ B 



P3 



2ir q 3 + q 3 _£ B 



4 



r 3 



' 3 d\ ' 
2tt q 3 + q 3 _£% 



r-i = 



2tt gH% 
3 d 3 



q+q- 



qi + q 2 - 



q\ qt d\ 



B 



(89) 
(90) 

(91) 



48' 
(92) 



,(r-r')+[qlv c (O)+q 2 v c (O)]/2 



(80) 



The coefficients in the fugacity 



/ e -*aq a <j>(r)-t/3qf,<l>(r')\ _ Q - [g* ti D H(0)+g£i>DH(0)] /2 
v p -^g a ?(9i*H(r-r') 



2 fr"lp- !Q 9»*( r ) 



<*V)e 



-9>dh(0)/2 



(81) 



«DH(0)-g>2 H ( r _ r ') 



(82) 



The brackets (• • •) denote averages over the fields 
and tp a where w _1 and v^, H are the propagators. 



In Eq. (|39|) we define the fugacity of negative ions A_ as an 

expansion in terms of A+. The coefficients ao, a%, etc. are 
determined by the condition of global electroneutrality. 
They are given by 

ao = — , (93) 



q+ r 2 2 1/ n ^B 



i + q - 

q+- 



(94) 



The coefficients in the grand-canonical free energy 

We give here the explicit expressions for the coefficients 
g, Eq. @. Using the function H(x), defined in Eq.@, 
the coefficients read 



TOl ^3^T^" 



In 



(83) 



a 2 = ■^-[-2m2q + +2mig_+m 3 [oi-g 2 ]]- 



6 q-d 3 . 



[q++q-Y 



x [7q 3 + -9g2 fl _ + 2 3] 



nqii 



2 «3 



! -[q\-q\q 2 -+q 2 +ql-q-] 



3 q-d\ 



xln(^g + [g+ +</-]), (95) 
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04 



24q[ 
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X {24m [5q+q 2 _ - q 3 _] + 24n 2 [q% - 5q+q 2 _] 

+72n 3 [q 2 + q--q+ql] +24mi [q\ql -iq+qt] +24m 2 [-4g* g_ 
- 4d + 7<£<£] + 12m 3 [2g 4 g_ _ q \ q 2 _ _ 6g 2 g3 + 5(?+(? 4 j 



7^ 



B [q+ +q-Y [26^ ~ 34, ?+9- ~ *W+<£ + 67^(7 



3 „2 1 fi7 „2 „3 



4 



- 45 W_ + 17g* ] } + ^ ^ [g + - ff _] [ S+ + g_] 5/2 

x [104 - llgij.g- + 13g^g 2 _ 4g+(? 3 + 3g 4 j 

/ 47rf 3 \ 
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3/0 5/2*9/2 + 
7T 3 / 2 0+' * fl ' r 1 r , n 5/2r in 4 



a 5 



6 q _dT 



\q + -q-][q + +q-Y' Z [lOq\-llqlq- 



+ 13q 2 + q 2 _ - 4q+q 3 _ + 3qi] In (-^9+ [q+ +«-])- (98) 
Notice that mi, m 2 , etc. were defined in Eqs. (|83|-pr 



The coefficients in the mean activity coefficient 

In Eq. [73] we obtained the low-density expansion of the 
mean activity coefficient of ionic solutions where the ions 
have valences q+ and g_ and effective diameters d+ and d_ 
(in Angstroms). Defining uj = 6.022 x 10~ 4 £pq^_q ( with 
77 = g_/g + and £ = d-/d+), the coefficients in Eq. |73| read 



1/1 = 



~KU>rj 
6[1 + »/] 



77^(1 + 7?], (99) 
{-l+3?7+87/ 2 +377 3 -7; 4 +4[i7(e + )+ln(e + ) 



-4?7 4 



— 877^ 



H - 



1 + C 



-hi 



.2776+ , 

1 + ^ 



+ 8? ? 2 [1 - ?7 2 ] \n(rj) - [2 - 4?7 2 + 2t] 4 ] ln(367ru;[l + 7?]) j 

7TCJ77 



^llog 



■[i-ri\ J [i + v], 



3/2 



|42— 677— 14t7 2 + 68t7 3 — I477 4 — 6?7 5 +4277 6 



24VT~R7 



4 0[1- f ] |ff(e+) +ln(e+)J + 4 77 6 [l - -] [h( , 



In 



II277 3 
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- 16?7 3 [7 - 2?7 2 + 5?7 3 ] In(r?) + 877 [l - rj 2 } 2 ln(367rw[l + 77]) 
-20[l + ?7 3 ] 2 ln(647rw[l + 77])|, (102) 
77 [710;] 3 / 2 , 

^3/21og = g V 1 + 7 / 

x [5 - 777 + 7?7 2 + 7t7 3 - 7?7 4 + 5?7 5 ] . (103) 
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